# KKT check the output theta matrix
kktchk <- function(obj, pf, thr) {
    lambda <- obj$lambda
    nk <- nrow(obj$delta)
    nvars <- ncol(obj$sigma_old)
    for (l in 1:length(lambda)) {
        cat("now checking lambda ", l, "\n")
        theta <- t(obj$theta[[l]])
        thetaInner <- apply(theta, 2, crossprod)
        thetaNorm <- sqrt(thetaInner)
        
        sdiag <- diag(obj$sigma_old)
        sjj <- t(replicate(nk, sdiag))
        thetaTilda <- (obj$delta - theta %*% obj$sigma_old + sjj * theta)/sjj
        thetaDif <- theta - thetaTilda
        
        for (j in 1:nvars) {
            los <- lambda[l] * pf[j]/obj$sigma_old[j, j]
            if (thetaNorm[j] == 0) {
                dif_norm <- sqrt(crossprod(thetaDif[, j]))
                tmp <- dif_norm - los
                if (tmp > 0) 
                  cat("violated at t > 0", tmp, "\n")
            } else {
                tmp3 <- thetaDif[, j] + los * theta[, j]/thetaNorm[j]
                if (any(abs(tmp3) > thr)) 
                  cat("violated at t = 0", tmp3, "\n")
            }
        }
    }
}

# extract fortran outputs and format it into sparse matries
formatoutput <- function(fit, maxit, pmax, nvars, vnames, nk) {
    nalam <- fit$nalam
    ntheta <- fit$ntheta[seq(nalam)]
    nthetamax <- max(ntheta)
    lam <- fit$alam[seq(nalam)]
    stepnames <- paste("s", seq(nalam) - 1, sep = "")
    resnames <- paste("delta", seq(nk), sep = "")
    
    errmsg <- err(fit$jerr, maxit, pmax)  ### error messages from fortran
    switch(paste(errmsg$n), `1` = stop(errmsg$msg, call. = FALSE), `-1` = cat(errmsg$msg))
    
    dd <- c(nvars, nk)
    df <- rep(0, nalam)
    if (nthetamax > 0) {
        ja <- fit$itheta[seq(nthetamax)]
        oja <- order(ja)
        ja <- rep(ja[oja], nk)
        itheta <- cumsum(c(1, rep(nthetamax, nk)))
        pos <- rep(1:nalam, each = nk * pmax)
        theta <- split(fit$theta[seq(nk * pmax * nalam)], pos)
        for (l in 1:nalam) {
            theta[[l]] <- matrix(theta[[l]], pmax, nk, byrow = TRUE)[seq(nthetamax), , drop = FALSE]
            df[l] <- sum(apply(theta[[l]], 1, sum) != 0)
            theta[[l]] <- new("dgCMatrix", Dim = dd, Dimnames = list(vnames, resnames), x = as.vector(theta[[l]][oja, 
                ]), p = as.integer(itheta - 1), i = as.integer(ja - 1))
        }
    } else {
        theta <- list()
        for (l in 1:nalam) {
            theta[[l]] <- zeromat(nvars, nk, vnames, resnames)
        }
        df <- rep(0, nalam)
    }
    list(theta = theta, df = df, dim = dd, lambda = lam)
}


# generate sigma, delta and mu from x, y
mlda.prep <- function(x, y) {
    # data setup
    x <- as.matrix(x)
    y <- drop(y)
    nclass <- as.integer(length(unique(y)))
	prior<-rep(0,nclass)
    for (k in 1:nclass) {
        prior[k] <- mean(y == k)
    }
    nvars <- as.integer(ncol(x))
    nobs <- nrow(x)
    nres <- length(y)
    if (nres != nobs) 
        stop("x and y have different number of observations")
    # prepare sigma and delta
    mu <- matrix(0, nvars, nclass)
    sigma <- matrix(0, nvars, nvars)
    for (i in 1:nclass) {
        mu[, i] <- apply(x[y == i, ], 2, mean)
        sigma <- sigma + (sum(y == i) - 1) * cov(x[y == i, ])
    }
    sigma <- sigma/(n - nclass)
    delta <- sweep(mu[, -1], 1, mu[, 1], "-")
    delta <- t(delta)
    outlist <- list(sigma = sigma, delta = delta, mu = mu, prior = prior)
    outlist
}

# convert theta list into theta array
list2array <- function(thetaList) {
    nlambda <- length(thetaList)
    nvars <- nrow(thetaList[[1]])
    nk <- ncol(thetaList[[1]])
    thetaArray <- array(0, dim = c(nk, nvars, nlambda))
    for (l in 1:nlambda) {
        thetaArray[, , l] <- t(as.matrix(thetaList[[l]]))
    }
    thetaArray
}


## These functions below are minor modifications or directly copied from the glmnet package:
## Jerome Friedman, Trevor Hastie, Robert Tibshirani (2010).  Regularization Paths for
## Generalized Linear Models via Coordinate Descent.  Journal of Statistical Software, 33(1),
## 1-22.  URL http://www.jstatsoft.org/v33/i01/.  The reason they are copied here is because
## they are internal functions and hence are not exported into the global environment.  The
## original comments and header are preserved.

err <- function(n, maxit, pmax) {
    if (n == 0) 
        msg <- ""
    if (n > 0) {
        # fatal error
        if (n < 7777) 
            msg <- "Memory allocation error; contact package maintainer"
        if (n == 10000) 
            msg <- "All penalty factors are <= 0"
        n <- 1
        msg <- paste("in the fortran code -", msg)
    }
    if (n < 0) {
        # non fatal error
        if (n > -10000) 
            msg <- paste("Convergence for ", -n, "th lambda value not reached after maxit=", maxit, 
                " iterations; solutions for larger lambdas returned", sep = "")
        if (n < -10000) 
            msg <- paste("Number of nonzero coefficients along the path exceeds pmax=", pmax, " at ", 
                -n - 10000, "th lambda value; solutions for larger lambdas returned", sep = "")
        n <- -1
        msg <- paste("from the fortran code -", msg)
    }
    list(n = n, msg = msg)
}

zeromat <- function(nvars, nalam, vnames, stepnames) {
    ca <- rep(0, nalam)
    ia <- seq(nalam + 1)
    ja <- rep(1, nalam)
    dd <- c(nvars, nalam)
    new("dgCMatrix", Dim = dd, Dimnames = list(vnames, stepnames), x = as.vector(ca), p = as.integer(ia - 
        1), i = as.integer(ja - 1))
}

lamfix <- function(lam) {
    llam <- log(lam)
    lam[1] <- exp(2 * llam[2] - llam[3])
    lam
} 
